#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

#include "merge.h"

//#define CUTOFF 1024

int rank(double x, double X[], int n) {
  int L = 0;
  int R = n;
  while (L < R) {
    int i = (L + R) / 2;
    if (X[i] < x) {
      L = i + 1;
    } else {
      R = i;
    }
  }
  return L;
}

void merge_divconq(double a[], long n, double b[], long m, double c[], long cutoff) {
  int i;

  //printf("%i %li %li %li\n", omp_get_thread_num(), cutoff, n, m);

  if (n == 0) { // task parallelize for large n
    for (i = 0; i < m; i++) c[i] = b[i];
  } else if (m == 0) { // task parallelize for large m
    for (i = 0; i < n; i++) c[i] = a[i];
  } else if (n + m <= cutoff) {
    //merge(a, n, b, m, c); // sequential merge for small problems
    seq_merge1(a, n, b, m, c);
  } else {
    int r = n / 2;
    int s = rank(a[r], b, m);
    c[r + s] = a[r];
    #pragma omp task
    merge_divconq(a, r, b, s, c, cutoff);
    #pragma omp task
    merge_divconq(&a[r + 1], n - r - 1, &b[s], m - s, &c[r + s + 1], cutoff);
  }
}

void merge(double a[], long n, double b[], long m, double c[]) {
  // replace this by a parallel merge algorithm
  //seq_merge1(a, n, b, m, c);

  #pragma omp parallel
  #pragma omp master
  {
    if (n >= m)
      merge_divconq(a, n, b, m, c, (n + m) / omp_get_max_threads());
    else
      merge_divconq(b, m, a, n, c, (n + m) / omp_get_max_threads());
  }
}
